library(data.table)
library(magrittr)
library(ggplot2)
theme_set(theme_minimal())

metadata <- fread("ihmp/data/ibdmdb_healthy.csv")
sra_data <- fread("ihmp/data/ibdmd_healthy_runtable.tsv")[, .(Run, `Library Name`)]
sra_data[, "External ID" := tstrsplit(`Library Name`, "_MGX")[[1]]]
metadata <- metadata[sra_data, on="External ID", nomatch=0]
metadata
man <- fread("db/data/manifest.csv")[, .(db, rank, matched_taxid, seqlength)]

foods <- fread("ihmp/data/food_abundance.csv")
foods <- man[foods, on="matched_taxid"]
foods[, "tpm" := 1.0e6 * reads / as.double(seqlength)]
foods[is.na(tpm), tpm := 0.0]

contents <- fread("ihmp/data/food_content.csv")
contents
tests <- list(
  vegetables = list(group="Vegetables (salad, tomatoes, onions, greens, carrots, peppers, green beans, etc)", ex=expression(food_group == "Vegetables")),
  fruits = list(group="Fruits (no juice) (Apples, raisins, bananas, oranges, strawberries, blueberries", ex=expression(food_group == "Fruits")),
  beans = list(group="Beans (tofu, soy, soy burgers, lentils, Mexican beans, lima beans etc)", ex=expression(food_subgroup == "Beans")),
  `white meat` = list(group="White meat (chicken, turkey, etc.)", ex=expression(food_subgroup == "Poultry")),
  shellfish = list(group="Shellfish (shrimp, lobster, scallops, etc.)", ex=expression(food_subgroup %chin% c("Mollusks", "Crustaceans"))),
  fish = list(group="Fish (fish nuggets, breaded fish, fish cakes, salmon, tuna, etc.)", ex=expression(food_subgroup == "Fishes")),
  `red meat` = list(group="Red meat (beef, hamburger, pork, lamb)", ex=expression(food_subgroup %chin% c("Swine", "Bovines", "Caprae")))
)
library(patchwork)

tables <- list()

freqs <- c(
  `No, I did not consume these products in the last 7 days` = 0,
  `Within the past 4 to 7 days` = 1/5.5,
  `Within the past 2 to 3 days` = 1/2.5,
  `Yesterday, 1 to 2 times` = 1.5,                                
  `Yesterday, 3 or more times` = 3,
  `NA` = 0
)

results <- list()
plots <- list()
merged <- metadata[foods, on=c(`Run` = "sample_id")]
for (i in 1:length(tests)) {
  name <- names(tests)[i]
  vals <- tests[[i]]
  dt <- merged[, .(abundance=sum(reads[eval(vals$ex)]), group=.SD[[vals$group]][1], total_reads=total_reads[1]), by="External ID"]
  ns <- dt[, .N, by="group"] |> setkey(group)
  dt <- dt[group %chin% names(freqs) & ns[group, N] >= 10]
  dt[, "frequency" := freqs[group]]
  dt[, group := factor(group, levels=names(freqs)[names(freqs) %in% group])]
  dt[, "relative" := (abundance + 1) / (total_reads + 1)]
  dt[, "item" := name]
  plots[[name]] <- ggplot(dt) + aes(x=group, y=log10(relative)) + 
    geom_jitter(width=0.3) + 
    geom_boxplot(width=0.1) + 
    stat_smooth(aes(group=1), method="lm") + 
    labs(title=name, x="consumption frequency [servings/day]", y="abundance")
  print(name)
  abmod <- lm(log10(relative) ~ frequency, data=dt)
  premod <- glm((abundance > 0) ~ frequency, data=dt)
  tables[[name]] <- dt
  results[[name]] <- dt[, .(n=.N, relative=mean(log10(relative*1e6)), sd=sd(log10(relative*1e6)), detected=sum(abundance > 0), item=name), by="group"]
  results[[name]][, "p_abundance" := summary(abmod)$coefficients[2, 4]]
  results[[name]][, "p_prevalence" := summary(premod)$coefficients[2, 4]]
}
[1] "vegetables"
[1] "fruits"
[1] "beans"
[1] "white meat"
[1] "shellfish"
[1] "fish"
[1] "red meat"
r <- rbindlist(results)
ta <- rbindlist(tables)
ggplot(r) +
  aes(x=group, y=relative) +
  geom_jitter(data=ta, aes(y=log10(relative*1e6)), stroke=0, size=0.5, width=0.3) +
  #geom_boxplot(outlier.color="NA") +
  geom_pointrange(aes(ymin=relative-sd, ymax=relative+sd), shape=21, fill="white") +
  geom_text(data=unique(r, by="item"), aes(x=0, y=Inf, label=sprintf("p[t-test]=%.2g\np[logit]=%.2g", p_abundance, p_prevalence)), size=3, vjust=1, hjust=0, nudge_y=0.5) +
  facet_wrap(~ item, ncol=2) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  labs(x="", y="log-relative abundance [log₁₀RPM]") +
  guides(color=F)

ggsave("figures/ihmp_foods.png", width=3, height=11, dpi=300)
rare <- foods[, .(n=sum(reads > 0), reads=sum(reads), total_reads=max(total_reads)), by="sample_id"]
ggplot(rare[n>0], aes(x=total_reads, y=reads)) + geom_point() + stat_smooth(method="lm", color="tomato") +
  labs(x="library size (#reads in sample)", y="#food reads") + scale_x_log10() + scale_y_log10()

cor.test(rare$total_reads, rare$reads)

    Pearson's product-moment correlation

data:  rare$total_reads and rare$reads
t = -0.16919, df = 363, p-value = 0.8657
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.11142887  0.09385626
sample estimates:
         cor 
-0.008879863 
ggsave("figures/food_reads.pdf", width=5, height=3)
ggplot(foods[seqlength>1e6], aes(x=seqlength/1e9, y=relative_abundance)) + geom_point() + stat_smooth(method="lm", color="tomato") +
  labs(x="haploid genome size [Gbps]", y="relative food abundance") + scale_x_log10() + scale_y_log10()

foods[seqlength>1e6, cor.test(reads, seqlength)]

    Pearson's product-moment correlation

data:  reads and seqlength
t = 1.9748, df = 2209, p-value = 0.04841
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.000294607 0.083521188
sample estimates:
       cor 
0.04198072 
ggsave("figures/genome_size.pdf", width=5, height=3)

Diet distances

species <- fread("ihmp/data/S_counts.csv")[grepl("Bacteria|Archaea", lineage)]
microbiome_mat <- dcast(species, sample_id ~ lineage, value.var = "new_est_reads", fill=0, fun.aggregate = sum)
sids <- microbiome_mat[, tstrsplit(sample_id, "S_")[[2]]]
microbiome_mat <- as.matrix(microbiome_mat[, "sample_id" := NULL])
rownames(microbiome_mat) <- sids
good <- colMeans(microbiome_mat) > 10
microbiome_mat <- microbiome_mat[, good]
# micro_rare <- phyloseq(otu_table(microbiome_mat, taxa_are_rows = F)) |> rarefy_even_depth(100000) |> otu_table()
microbiome_relative <- microbiome_mat / rowSums(microbiome_mat)
diet <- metadata[, names(metadata)[72:92][-11], with=F]
diet <- apply(diet, 2, function(x) ifelse(x == "", 0, freqs[x]))
rownames(diet) <- metadata[["Run"]]
diet <- diet[rowSums(diet) > 0, ]
diet <- diet[, colMeans(diet >0) > 0.1]
dim(diet)
[1] 362  19
foodab <- foods
foodab[, "id" := paste(matched_taxid, species, wikipedia_id, sep="|")]
food_mat <- dcast(foodab, sample_id ~ id, value.var = "reads", fill=0, fun.aggregate = sum)
sids <- food_mat[, sample_id]
food_mat <- as.matrix(food_mat[, sample_id := NULL])
rownames(food_mat) <- sids
food_relative <- food_mat[rowSums(food_mat) > 0, ]
food_relative <- food_relative / rowSums(food_relative)
food_relative <- food_relative[rowSums(food_relative) > 0, ]
food_relative <- food_relative[, colMeans(food_relative >0) > 0.1]
dim(food_relative)
[1] 359  16

Beta diversity tests

Let’s write a little function that runs the test for microbiome <-> other comparison and plots and returns the results.

library(vegan)
Loading required package: permute
Loading required package: lattice
This is vegan 2.6-8
micro_mantel <- function(first, other, description) {
  sids <- intersect(rownames(first), rownames(other))
  micro_dist <- vegdist(first[sids, ], "bray")
  other_dist <- vegdist(other[sids, ], "bray")
  
  test <- mantel(micro_dist, other_dist, method="pearson")
  results <- data.table(micro=as.numeric(micro_dist), other=as.numeric(other_dist))
  pl <- ggplot(results) +
    aes(x=other, y=micro) +
    geom_point(size=1, alpha=0.25, stroke=0) +
    stat_smooth(method="lm") +
    labs(x=paste(description, "[Bray]"), y="species abundance distance [Bray]")
  print(pl)
  print(test)
}

Now we run it for the measures.

micro_mantel(microbiome_relative, diet, "FFQ distances")

Mantel statistic based on Pearson's product-moment correlation 

Call:
mantel(xdis = micro_dist, ydis = other_dist, method = "pearson") 

Mantel statistic r: 0.04132 
      Significance: 0.079 

Upper quantiles of permutations (null model):
   90%    95%  97.5%    99% 
0.0356 0.0482 0.0558 0.0649 
Permutation: free
Number of permutations: 999

micro_mantel(microbiome_relative, food_relative, "MEDI food distances")
Warning: you have empty rows: their dissimilarities may be
                 meaningless in method “bray”

Mantel statistic based on Pearson's product-moment correlation 

Call:
mantel(xdis = micro_dist, ydis = other_dist, method = "pearson") 

Mantel statistic r: 0.08119 
      Significance: 0.001 

Upper quantiles of permutations (null model):
   90%    95%  97.5%    99% 
0.0273 0.0335 0.0399 0.0461 
Permutation: free
Number of permutations: 999

micro_mantel(food_relative, diet, "FFQ vs. MEDI")
Warning: you have empty rows: their dissimilarities may be
                 meaningless in method “bray”

Mantel statistic based on Pearson's product-moment correlation 

Call:
mantel(xdis = micro_dist, ydis = other_dist, method = "pearson") 

Mantel statistic r: 0.08161 
      Significance: 0.002 

Upper quantiles of permutations (null model):
   90%    95%  97.5%    99% 
0.0307 0.0409 0.0476 0.0585 
Permutation: free
Number of permutations: 999

LS0tCnRpdGxlOiAiaUhNUCIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKYGBge3J9CmxpYnJhcnkoZGF0YS50YWJsZSkKbGlicmFyeShtYWdyaXR0cikKbGlicmFyeShnZ3Bsb3QyKQp0aGVtZV9zZXQodGhlbWVfbWluaW1hbCgpKQoKbWV0YWRhdGEgPC0gZnJlYWQoImlobXAvZGF0YS9pYmRtZGJfaGVhbHRoeS5jc3YiKQpzcmFfZGF0YSA8LSBmcmVhZCgiaWhtcC9kYXRhL2liZG1kX2hlYWx0aHlfcnVudGFibGUudHN2IilbLCAuKFJ1biwgYExpYnJhcnkgTmFtZWApXQpzcmFfZGF0YVssICJFeHRlcm5hbCBJRCIgOj0gdHN0cnNwbGl0KGBMaWJyYXJ5IE5hbWVgLCAiX01HWCIpW1sxXV1dCm1ldGFkYXRhIDwtIG1ldGFkYXRhW3NyYV9kYXRhLCBvbj0iRXh0ZXJuYWwgSUQiLCBub21hdGNoPTBdCm1ldGFkYXRhCmBgYAoKYGBge3J9Cm1hbiA8LSBmcmVhZCgiZGIvZGF0YS9tYW5pZmVzdC5jc3YiKVssIC4oZGIsIHJhbmssIG1hdGNoZWRfdGF4aWQsIHNlcWxlbmd0aCldCgpmb29kcyA8LSBmcmVhZCgiaWhtcC9kYXRhL2Zvb2RfYWJ1bmRhbmNlLmNzdiIpCmZvb2RzIDwtIG1hbltmb29kcywgb249Im1hdGNoZWRfdGF4aWQiXQpmb29kc1ssICJ0cG0iIDo9IDEuMGU2ICogcmVhZHMgLyBhcy5kb3VibGUoc2VxbGVuZ3RoKV0KZm9vZHNbaXMubmEodHBtKSwgdHBtIDo9IDAuMF0KCmNvbnRlbnRzIDwtIGZyZWFkKCJpaG1wL2RhdGEvZm9vZF9jb250ZW50LmNzdiIpCmNvbnRlbnRzCmBgYAoKYGBge3J9CnRlc3RzIDwtIGxpc3QoCiAgdmVnZXRhYmxlcyA9IGxpc3QoZ3JvdXA9IlZlZ2V0YWJsZXMgKHNhbGFkLCB0b21hdG9lcywgb25pb25zLCBncmVlbnMsIGNhcnJvdHMsIHBlcHBlcnMsIGdyZWVuIGJlYW5zLCBldGMpIiwgZXg9ZXhwcmVzc2lvbihmb29kX2dyb3VwID09ICJWZWdldGFibGVzIikpLAogIGZydWl0cyA9IGxpc3QoZ3JvdXA9IkZydWl0cyAobm8ganVpY2UpIChBcHBsZXMsIHJhaXNpbnMsIGJhbmFuYXMsIG9yYW5nZXMsIHN0cmF3YmVycmllcywgYmx1ZWJlcnJpZXMiLCBleD1leHByZXNzaW9uKGZvb2RfZ3JvdXAgPT0gIkZydWl0cyIpKSwKICBiZWFucyA9IGxpc3QoZ3JvdXA9IkJlYW5zICh0b2Z1LCBzb3ksIHNveSBidXJnZXJzLCBsZW50aWxzLCBNZXhpY2FuIGJlYW5zLCBsaW1hIGJlYW5zIGV0YykiLCBleD1leHByZXNzaW9uKGZvb2Rfc3ViZ3JvdXAgPT0gIkJlYW5zIikpLAogIGB3aGl0ZSBtZWF0YCA9IGxpc3QoZ3JvdXA9IldoaXRlIG1lYXQgKGNoaWNrZW4sIHR1cmtleSwgZXRjLikiLCBleD1leHByZXNzaW9uKGZvb2Rfc3ViZ3JvdXAgPT0gIlBvdWx0cnkiKSksCiAgc2hlbGxmaXNoID0gbGlzdChncm91cD0iU2hlbGxmaXNoIChzaHJpbXAsIGxvYnN0ZXIsIHNjYWxsb3BzLCBldGMuKSIsIGV4PWV4cHJlc3Npb24oZm9vZF9zdWJncm91cCAlY2hpbiUgYygiTW9sbHVza3MiLCAiQ3J1c3RhY2VhbnMiKSkpLAogIGZpc2ggPSBsaXN0KGdyb3VwPSJGaXNoIChmaXNoIG51Z2dldHMsIGJyZWFkZWQgZmlzaCwgZmlzaCBjYWtlcywgc2FsbW9uLCB0dW5hLCBldGMuKSIsIGV4PWV4cHJlc3Npb24oZm9vZF9zdWJncm91cCA9PSAiRmlzaGVzIikpLAogIGByZWQgbWVhdGAgPSBsaXN0KGdyb3VwPSJSZWQgbWVhdCAoYmVlZiwgaGFtYnVyZ2VyLCBwb3JrLCBsYW1iKSIsIGV4PWV4cHJlc3Npb24oZm9vZF9zdWJncm91cCAlY2hpbiUgYygiU3dpbmUiLCAiQm92aW5lcyIsICJDYXByYWUiKSkpCikKYGBgCgoKCmBgYHtyLCBmaWcud2lkdGg9NCwgZmlnLmhlaWdodD0zfQpsaWJyYXJ5KHBhdGNod29yaykKCnRhYmxlcyA8LSBsaXN0KCkKCmZyZXFzIDwtIGMoCiAgYE5vLCBJIGRpZCBub3QgY29uc3VtZSB0aGVzZSBwcm9kdWN0cyBpbiB0aGUgbGFzdCA3IGRheXNgID0gMCwKICBgV2l0aGluIHRoZSBwYXN0IDQgdG8gNyBkYXlzYCA9IDEvNS41LAogIGBXaXRoaW4gdGhlIHBhc3QgMiB0byAzIGRheXNgID0gMS8yLjUsCiAgYFllc3RlcmRheSwgMSB0byAyIHRpbWVzYCA9IDEuNSwgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIAogIGBZZXN0ZXJkYXksIDMgb3IgbW9yZSB0aW1lc2AgPSAzLAogIGBOQWAgPSAwCikKCnJlc3VsdHMgPC0gbGlzdCgpCnBsb3RzIDwtIGxpc3QoKQptZXJnZWQgPC0gbWV0YWRhdGFbZm9vZHMsIG9uPWMoYFJ1bmAgPSAic2FtcGxlX2lkIildCmZvciAoaSBpbiAxOmxlbmd0aCh0ZXN0cykpIHsKICBuYW1lIDwtIG5hbWVzKHRlc3RzKVtpXQogIHZhbHMgPC0gdGVzdHNbW2ldXQogIGR0IDwtIG1lcmdlZFssIC4oYWJ1bmRhbmNlPXN1bShyZWFkc1tldmFsKHZhbHMkZXgpXSksIGdyb3VwPS5TRFtbdmFscyRncm91cF1dWzFdLCB0b3RhbF9yZWFkcz10b3RhbF9yZWFkc1sxXSksIGJ5PSJFeHRlcm5hbCBJRCJdCiAgbnMgPC0gZHRbLCAuTiwgYnk9Imdyb3VwIl0gfD4gc2V0a2V5KGdyb3VwKQogIGR0IDwtIGR0W2dyb3VwICVjaGluJSBuYW1lcyhmcmVxcykgJiBuc1tncm91cCwgTl0gPj0gMTBdCiAgZHRbLCAiZnJlcXVlbmN5IiA6PSBmcmVxc1tncm91cF1dCiAgZHRbLCBncm91cCA6PSBmYWN0b3IoZ3JvdXAsIGxldmVscz1uYW1lcyhmcmVxcylbbmFtZXMoZnJlcXMpICVpbiUgZ3JvdXBdKV0KICBkdFssICJyZWxhdGl2ZSIgOj0gKGFidW5kYW5jZSArIDEpIC8gKHRvdGFsX3JlYWRzICsgMSldCiAgZHRbLCAiaXRlbSIgOj0gbmFtZV0KICBwbG90c1tbbmFtZV1dIDwtIGdncGxvdChkdCkgKyBhZXMoeD1ncm91cCwgeT1sb2cxMChyZWxhdGl2ZSkpICsgCiAgICBnZW9tX2ppdHRlcih3aWR0aD0wLjMpICsgCiAgICBnZW9tX2JveHBsb3Qod2lkdGg9MC4xKSArIAogICAgc3RhdF9zbW9vdGgoYWVzKGdyb3VwPTEpLCBtZXRob2Q9ImxtIikgKyAKICAgIGxhYnModGl0bGU9bmFtZSwgeD0iY29uc3VtcHRpb24gZnJlcXVlbmN5IFtzZXJ2aW5ncy9kYXldIiwgeT0iYWJ1bmRhbmNlIikKICBwcmludChuYW1lKQogIGFibW9kIDwtIGxtKGxvZzEwKHJlbGF0aXZlKSB+IGZyZXF1ZW5jeSwgZGF0YT1kdCkKICBwcmVtb2QgPC0gZ2xtKChhYnVuZGFuY2UgPiAwKSB+IGZyZXF1ZW5jeSwgZGF0YT1kdCkKICB0YWJsZXNbW25hbWVdXSA8LSBkdAogIHJlc3VsdHNbW25hbWVdXSA8LSBkdFssIC4obj0uTiwgcmVsYXRpdmU9bWVhbihsb2cxMChyZWxhdGl2ZSoxZTYpKSwgc2Q9c2QobG9nMTAocmVsYXRpdmUqMWU2KSksIGRldGVjdGVkPXN1bShhYnVuZGFuY2UgPiAwKSwgaXRlbT1uYW1lKSwgYnk9Imdyb3VwIl0KICByZXN1bHRzW1tuYW1lXV1bLCAicF9hYnVuZGFuY2UiIDo9IHN1bW1hcnkoYWJtb2QpJGNvZWZmaWNpZW50c1syLCA0XV0KICByZXN1bHRzW1tuYW1lXV1bLCAicF9wcmV2YWxlbmNlIiA6PSBzdW1tYXJ5KHByZW1vZCkkY29lZmZpY2llbnRzWzIsIDRdXQp9CmBgYAoKYGBge3IsIGZpZy53aWR0aD0zLCBmaWcuaGVpZ2h0PTExfQpyIDwtIHJiaW5kbGlzdChyZXN1bHRzKQp0YSA8LSByYmluZGxpc3QodGFibGVzKQpnZ3Bsb3QocikgKwogIGFlcyh4PWdyb3VwLCB5PXJlbGF0aXZlKSArCiAgZ2VvbV9qaXR0ZXIoZGF0YT10YSwgYWVzKHk9bG9nMTAocmVsYXRpdmUqMWU2KSksIHN0cm9rZT0wLCBzaXplPTAuNSwgd2lkdGg9MC4zKSArCiAgI2dlb21fYm94cGxvdChvdXRsaWVyLmNvbG9yPSJOQSIpICsKICBnZW9tX3BvaW50cmFuZ2UoYWVzKHltaW49cmVsYXRpdmUtc2QsIHltYXg9cmVsYXRpdmUrc2QpLCBzaGFwZT0yMSwgZmlsbD0id2hpdGUiKSArCiAgZ2VvbV90ZXh0KGRhdGE9dW5pcXVlKHIsIGJ5PSJpdGVtIiksIGFlcyh4PTAsIHk9SW5mLCBsYWJlbD1zcHJpbnRmKCJwW3QtdGVzdF09JS4yZ1xucFtsb2dpdF09JS4yZyIsIHBfYWJ1bmRhbmNlLCBwX3ByZXZhbGVuY2UpKSwgc2l6ZT0zLCB2anVzdD0xLCBoanVzdD0wLCBudWRnZV95PTAuNSkgKwogIGZhY2V0X3dyYXAofiBpdGVtLCBuY29sPTIpICsKICB0aGVtZShheGlzLnRleHQueCA9IGVsZW1lbnRfdGV4dChhbmdsZSA9IDkwLCB2anVzdCA9IDAuNSwgaGp1c3Q9MSkpICsKICBsYWJzKHg9IiIsIHk9ImxvZy1yZWxhdGl2ZSBhYnVuZGFuY2UgW2xvZ+KCgeKCgFJQTV0iKSArCiAgZ3VpZGVzKGNvbG9yPUYpCmdnc2F2ZSgiZmlndXJlcy9paG1wX2Zvb2RzLnBuZyIsIHdpZHRoPTMsIGhlaWdodD0xMSwgZHBpPTMwMCkKYGBgCgoKYGBge3IsIGZpZy53aWR0aD01LCBmaWcuaGVpZ2h0PTN9CnJhcmUgPC0gZm9vZHNbLCAuKG49c3VtKHJlYWRzID4gMCksIHJlYWRzPXN1bShyZWFkcyksIHRvdGFsX3JlYWRzPW1heCh0b3RhbF9yZWFkcykpLCBieT0ic2FtcGxlX2lkIl0KZ2dwbG90KHJhcmVbbj4wXSwgYWVzKHg9dG90YWxfcmVhZHMsIHk9cmVhZHMpKSArIGdlb21fcG9pbnQoKSArIHN0YXRfc21vb3RoKG1ldGhvZD0ibG0iLCBjb2xvcj0idG9tYXRvIikgKwogIGxhYnMoeD0ibGlicmFyeSBzaXplICgjcmVhZHMgaW4gc2FtcGxlKSIsIHk9IiNmb29kIHJlYWRzIikgKyBzY2FsZV94X2xvZzEwKCkgKyBzY2FsZV95X2xvZzEwKCkKY29yLnRlc3QocmFyZSR0b3RhbF9yZWFkcywgcmFyZSRyZWFkcykKZ2dzYXZlKCJmaWd1cmVzL2Zvb2RfcmVhZHMucGRmIiwgd2lkdGg9NSwgaGVpZ2h0PTMpCmBgYCAKCmBgYHtyLCBmaWcud2lkdGg9NSwgZmlnLmhlaWdodD0zfQpnZ3Bsb3QoZm9vZHNbc2VxbGVuZ3RoPjFlNl0sIGFlcyh4PXNlcWxlbmd0aC8xZTksIHk9cmVsYXRpdmVfYWJ1bmRhbmNlKSkgKyBnZW9tX3BvaW50KCkgKyBzdGF0X3Ntb290aChtZXRob2Q9ImxtIiwgY29sb3I9InRvbWF0byIpICsKICBsYWJzKHg9ImhhcGxvaWQgZ2Vub21lIHNpemUgW0dicHNdIiwgeT0icmVsYXRpdmUgZm9vZCBhYnVuZGFuY2UiKSArIHNjYWxlX3hfbG9nMTAoKSArIHNjYWxlX3lfbG9nMTAoKQpmb29kc1tzZXFsZW5ndGg+MWU2LCBjb3IudGVzdChyZWFkcywgc2VxbGVuZ3RoKV0KZ2dzYXZlKCJmaWd1cmVzL2dlbm9tZV9zaXplLnBkZiIsIHdpZHRoPTUsIGhlaWdodD0zKQpgYGAgCgojIyBEaWV0IGRpc3RhbmNlcwoKYGBge3J9CnNwZWNpZXMgPC0gZnJlYWQoImlobXAvZGF0YS9TX2NvdW50cy5jc3YiKVtncmVwbCgiQmFjdGVyaWF8QXJjaGFlYSIsIGxpbmVhZ2UpXQptaWNyb2Jpb21lX21hdCA8LSBkY2FzdChzcGVjaWVzLCBzYW1wbGVfaWQgfiBsaW5lYWdlLCB2YWx1ZS52YXIgPSAibmV3X2VzdF9yZWFkcyIsIGZpbGw9MCwgZnVuLmFnZ3JlZ2F0ZSA9IHN1bSkKc2lkcyA8LSBtaWNyb2Jpb21lX21hdFssIHRzdHJzcGxpdChzYW1wbGVfaWQsICJTXyIpW1syXV1dCm1pY3JvYmlvbWVfbWF0IDwtIGFzLm1hdHJpeChtaWNyb2Jpb21lX21hdFssICJzYW1wbGVfaWQiIDo9IE5VTExdKQpyb3duYW1lcyhtaWNyb2Jpb21lX21hdCkgPC0gc2lkcwpnb29kIDwtIGNvbE1lYW5zKG1pY3JvYmlvbWVfbWF0KSA+IDEwCm1pY3JvYmlvbWVfbWF0IDwtIG1pY3JvYmlvbWVfbWF0WywgZ29vZF0KIyBtaWNyb19yYXJlIDwtIHBoeWxvc2VxKG90dV90YWJsZShtaWNyb2Jpb21lX21hdCwgdGF4YV9hcmVfcm93cyA9IEYpKSB8PiByYXJlZnlfZXZlbl9kZXB0aCgxMDAwMDApIHw+IG90dV90YWJsZSgpCm1pY3JvYmlvbWVfcmVsYXRpdmUgPC0gbWljcm9iaW9tZV9tYXQgLyByb3dTdW1zKG1pY3JvYmlvbWVfbWF0KQpgYGAKCmBgYHtyfQpkaWV0IDwtIG1ldGFkYXRhWywgbmFtZXMobWV0YWRhdGEpWzcyOjkyXVstMTFdLCB3aXRoPUZdCmRpZXQgPC0gYXBwbHkoZGlldCwgMiwgZnVuY3Rpb24oeCkgaWZlbHNlKHggPT0gIiIsIDAsIGZyZXFzW3hdKSkKcm93bmFtZXMoZGlldCkgPC0gbWV0YWRhdGFbWyJSdW4iXV0KZGlldCA8LSBkaWV0W3Jvd1N1bXMoZGlldCkgPiAwLCBdCmRpZXQgPC0gZGlldFssIGNvbE1lYW5zKGRpZXQgPjApID4gMC4xXQpkaW0oZGlldCkKYGBgCgpgYGB7cn0KZm9vZGFiIDwtIGZvb2RzCmZvb2RhYlssICJpZCIgOj0gcGFzdGUobWF0Y2hlZF90YXhpZCwgc3BlY2llcywgd2lraXBlZGlhX2lkLCBzZXA9InwiKV0KZm9vZF9tYXQgPC0gZGNhc3QoZm9vZGFiLCBzYW1wbGVfaWQgfiBpZCwgdmFsdWUudmFyID0gInJlYWRzIiwgZmlsbD0wLCBmdW4uYWdncmVnYXRlID0gc3VtKQpzaWRzIDwtIGZvb2RfbWF0Wywgc2FtcGxlX2lkXQpmb29kX21hdCA8LSBhcy5tYXRyaXgoZm9vZF9tYXRbLCBzYW1wbGVfaWQgOj0gTlVMTF0pCnJvd25hbWVzKGZvb2RfbWF0KSA8LSBzaWRzCmZvb2RfcmVsYXRpdmUgPC0gZm9vZF9tYXRbcm93U3Vtcyhmb29kX21hdCkgPiAwLCBdCmZvb2RfcmVsYXRpdmUgPC0gZm9vZF9yZWxhdGl2ZSAvIHJvd1N1bXMoZm9vZF9yZWxhdGl2ZSkKZm9vZF9yZWxhdGl2ZSA8LSBmb29kX3JlbGF0aXZlW3Jvd1N1bXMoZm9vZF9yZWxhdGl2ZSkgPiAwLCBdCmZvb2RfcmVsYXRpdmUgPC0gZm9vZF9yZWxhdGl2ZVssIGNvbE1lYW5zKGZvb2RfcmVsYXRpdmUgPjApID4gMC4xXQpkaW0oZm9vZF9yZWxhdGl2ZSkKYGBgCgojIyBCZXRhIGRpdmVyc2l0eSB0ZXN0cwoKTGV0J3Mgd3JpdGUgYSBsaXR0bGUgZnVuY3Rpb24gdGhhdCBydW5zIHRoZSB0ZXN0IGZvciBtaWNyb2Jpb21lIDwtPiBvdGhlciBjb21wYXJpc29uIGFuZApwbG90cyBhbmQgcmV0dXJucyB0aGUgcmVzdWx0cy4KCmBgYHtyfQpsaWJyYXJ5KHZlZ2FuKQoKbWljcm9fbWFudGVsIDwtIGZ1bmN0aW9uKGZpcnN0LCBvdGhlciwgZGVzY3JpcHRpb24pIHsKICBzaWRzIDwtIGludGVyc2VjdChyb3duYW1lcyhmaXJzdCksIHJvd25hbWVzKG90aGVyKSkKICBtaWNyb19kaXN0IDwtIHZlZ2Rpc3QoZmlyc3Rbc2lkcywgXSwgImJyYXkiKQogIG90aGVyX2Rpc3QgPC0gdmVnZGlzdChvdGhlcltzaWRzLCBdLCAiYnJheSIpCiAgCiAgdGVzdCA8LSBtYW50ZWwobWljcm9fZGlzdCwgb3RoZXJfZGlzdCwgbWV0aG9kPSJwZWFyc29uIikKICByZXN1bHRzIDwtIGRhdGEudGFibGUobWljcm89YXMubnVtZXJpYyhtaWNyb19kaXN0KSwgb3RoZXI9YXMubnVtZXJpYyhvdGhlcl9kaXN0KSkKICBwbCA8LSBnZ3Bsb3QocmVzdWx0cykgKwogICAgYWVzKHg9b3RoZXIsIHk9bWljcm8pICsKICAgIGdlb21fcG9pbnQoc2l6ZT0xLCBhbHBoYT0wLjI1LCBzdHJva2U9MCkgKwogICAgc3RhdF9zbW9vdGgobWV0aG9kPSJsbSIpICsKICAgIGxhYnMoeD1wYXN0ZShkZXNjcmlwdGlvbiwgIltCcmF5XSIpLCB5PSJzcGVjaWVzIGFidW5kYW5jZSBkaXN0YW5jZSBbQnJheV0iKQogIHByaW50KHBsKQogIHByaW50KHRlc3QpCn0KYGBgCgpOb3cgd2UgcnVuIGl0IGZvciB0aGUgbWVhc3VyZXMuCgpgYGB7cn0KbWljcm9fbWFudGVsKG1pY3JvYmlvbWVfcmVsYXRpdmUsIGRpZXQsICJGRlEgdnMuIE1pY3JvYmlvbWUiKQpgYGAKCmBgYHtyfQptaWNyb19tYW50ZWwobWljcm9iaW9tZV9yZWxhdGl2ZSwgZm9vZF9yZWxhdGl2ZSwgIk1FREkgdnMuIE1pY3JvYmlvbWUiKQpgYGAKCmBgYHtyfQptaWNyb19tYW50ZWwoZm9vZF9yZWxhdGl2ZSwgZGlldCwgIkZGUSB2cy4gTUVESSIpCmBgYAo=